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Abstract 

We study three models of driven sandpile-type automata in the presence of 
quenched random defects. When the dynamics is conservative, all these models, 
termed the random sites (A), random bonds (B), and random slopes (C), self- 
organize into a critical state. For Model C the concentration-dependent exponents 
I/"") | are nonuniversal. In the case of nonconservative defects, the asymptotic state 

is subcritical. Possible defect-mediated nonequilibrium phase transitions are also 
"j^ . discussed. 

B 

1 Introduction 

o 
o 



Self-organized criticality (SOC) |lj is a dynamic phenomenon which occurs in certain 
dissipative systems with large numbers of degrees of freedom. Such a system, when slowly 
driven into its metastable state, self-organizes in a state with long-range correlations, 
similar to the critical state at a second-order phase transition. 

In conventional criticality, quenched disorder can be a relevant perturbation in the vicin- 
ity of the critical point. It is therefore natural to ask how the SOC state responds to 
similar perturbations. One might expect that self-organizing systems are more robust 
against random perturbations, since the critical state is an attractor of the dynamics, 
although the universality class may change in presence of disorder. 
We explore this question in the present work, where we demonstrate, using numerical 
simulations on simple models of self-organizing cellular automata with frozen random 
defects, the conditions for a system to self-organize, and determine the universality class 
of the critical behavior. Disorder-mediated phase transition between different types of 
metastable states are also discussed. 
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We study three kinds of random defects which locally affect the rules of relaxation in a 
manner analogous to the random site, random bond, and random field defects in spin 
models displaying conventional critical phenomena. All these models are based on the 
directed abelian 2-dimensional critical height model (toppling if h(i,j) > h c ) which is 
exactly solvable in the absence of defects ||, in which the dynamic rules are locally 
modified at a fraction c of "defect" sites. In model A, holes of infinite depth are placed 
at random sites. In model B there are two variables h\(i,j) and h 2 (i,j) associated with 
each site These are coupled at random sites (details are given in Section 3 below), 

in a manner so as to lead to a multiplicity of metastable states, as in the case of spin- 
glasses and other frustrated systems. In model C, frozen-in local slopes are introduced by 
preventing relaxation through the height instability at defect sites, and instead applying 
a critical slope toppling rule at all sites with slopes cr(i,j) > a c . In this case, regions with 
local slopes cr(i,j) ~ o c are established, while the rest of the system relaxes according to 
the critical height rule (the critical height h c is chosen such that h c <cr c ). In all three cases 
the ratio between the number of particles leaving one site and the number appearing at 
its neighbors is not uniform at defect sites. 

The numerical results presented here are from simulations on lattices (with periodic 
boundary conditions in the transverse direction) of size 12 < L < 384 for time steps (see 
below) up to 1 x 10 6 . A lattice with frozen-in defects is prepared and kept fixed for the 
entire number of time steps, and then a new configuration is prepared; results are then 
averaged over the total number of configurations. In each case, the system was driven 
by randomly adding particles to sites at the top (first row), h(l,j) — > h(l,j) + 1. 

2 Subcriticality: random-site model (A) 

The dynamic height variable h(i,j), associated with each lattice site on a square 
lattice, is updated according to the rules of the critical height model: if h(i,j) exceeds a 
critical value h c , then the site is unstable and relaxes according to 

h(i,j) -+h{i,j)-2 ; h(i + l,j±)^h(i + l,j±) + l, (1) 

where (z + l,j±) are the two neighboring downstream sites. This rule applies to all sites 
except for a fraction c of randomly distributed defect sites, at which infinitely deep holes 
are placed. Thus two grains are lost each time an avalanche hits a defect, rendering the 
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dynamic process nonconservative. (Similar nonconservative models have been considered 
earlier ||). For c = the critical exponents are listed in Table 1. An annealed version of 
defects of the above type was studied numerically in Ref. where it was shown that the 
system self-organizes into a subcritical state. The probability distributions of duration 
P(t, c) and size D(s, c) were found to fulfill the following scaling forms 

Pit, c) = r"V{t/Z t (c)) ; D(s, c) = s- T V(s/Uc)) , (2) 

where the correlation lengths £t(c) ~ c~ M * and £ s (c) ~ c _Ms for time and spatial corre- 
lations, respectively, diverge in the limit c — > with fi t = 1 and fi s = 1.5 ||, |]. The 
scaling functions defined in Eq. (H) were calculated numerically in Ref. M. With the 
dynamical rules (jj) there are no recurrences in the avalanche dynamics and thus both 
annealed and quenched (frozen) impurities lead to the same universal scaling properties. 
In Fig. 1 we plot < n(l) >, the average number of particles relaxed || up to distance I 
(measured from the top of the pile) for various values of the defect concentration: 
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Figure 1: Time average of the number of topplings < n(l) > vs. distance I from the top 
of pile for various concentrations of defects c. 

In the absence of defects < n(l) > is proportional to the average number of toppled sites 
< s(l) > and thus exhibits a power-law behavior 



< n(l) >~ l» 



(3) 



with /i = 1 (cf. the top curve in Fig. 1). For finite concentration of defects c, the curve 
flattens, and the scaling region with slope fi decreases with increasing c, corresponding 
to the decreasing correlation length in the system. At a critical concentration c = c* = 
0.295 = 1 — pa, the directed-percolation threshold |7[], it is no longer possible to have a 
lattice spanning cluster M. 

For sandpiles without defects, the fluxes into and out of the lattice are equal- the dis- 
tribution G(J, L) of the current which flows over the rim of the system is shown in Fig. 
2 for c = 0.0 and various systems sizes L. As expected for a self-organized critical state, 
the outflow current distribution exhibits scale invariance [H, i.e. 



G(J,L) = L-Pg(J/L^) 



(4) 
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Figure 2: Outflow current distribution in the absence of defects for system sizes L = 24, 
48, 96, 192, and 384 (left), and the corresponding finite- size scaling plot (right). 



The finite-size scaling fit, also shown in Fig. 2 has the exponents (3 = 1 and (f) = 0.5. 
With defects, the outflow current diminishes with increasing concentration of defects, 
eventually vanishing at c = c* (see Section 5). There is no apparent scaling form of 
G(J,c). 
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3 Universal criticality: random-bond model (B) 



In order to simulate effects of random bonds in a sandpile automaton, we introduce a two- 
state variable h = (hi, h 2 ) at each lattice site where h\ and h 2 are not necessarily 
integer. Each bond carries a quenched variable b with value b = ±1: the disorder here 
is that a random fraction c of the bonds have b = — 1. The evolution rule for model B 
depends on the absolute value of the difference between components hi and h 2 at a site, 
which causes instability if it exceeds a critical value d c . The entire number of particles 
then topple, and the two downstream neighboring sites are updated as follows (similar 
models were discussed earlier in Ref. 0): 

If \h(i,j)-h 2 (i,j)\>d c (5) 

then 

h{i,j)^0 ; h 2 (i,j)^0 (6) 

and 

h 1 (i + l,j±)^h 1 (i + l,j±) + (\h 1 (i,j) + h 2 (i,j))/2, if b = +l (7) 

hv(i + l,j±)^h 2 (i + l,j±) + (Xh 2 (i,j) + h 1 (i,j))/2, if b = -l (8) 

For A = 1 the dynamics is conservative: all particles which leave site appear at 
its neighbors. The instability condition Eq. (J|) leads to a variety of states with high 
local values of hi and h 2 . The asymptotic state is, however, SOC, and it appears that 
the precise form of the coupling between two states hi and h 2 is unimportant (see also 
Ref. (l^l). Such a model incorporates some features of neural networks, and introduces 
frustration effects 0. 

It is necessary to have conservative dynamics (A = 1) in order for the system to self- 
organize into a critical state. Results of numerical simulations for the distributions of 
duration P(t) and size D(s > sq) are shown in Fig. 3 for concentration c = 0.5, A = 1 
(solid lines) and A = 0.9 (open circles). The present results average over simulations 
for a long time (~ 10 6 steps) from several configurations, each of which is kept fixed 
for a particular simulation. The slopes of the straight lines in the case of conservative 
dynamics (A = 1) determine the critical exponents according to Pit) ~ t~( 1+e ) and 
D(s > s ) ~ s~ T . We find the best fit for 6 = 1.040 ± 0.018 and r = 0.650 ± 0.028. 
In this model the number of relaxed "particles" n is not proportional to the number of 
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Figure 3: Distributions of duration P(t) and size D(s > sq) of the relaxation clusters 
in model B for conserving A = 1 (solid lines) and nonconserving dynamics A<1 (open 
circles) . 

sites s at which relaxation occurs, thus leading to a new distribution Q(n) ~ n~ t < l+Tn \ 
We also calculated <n(l)>, to obtain the exponent \i defined in Eq. (§), as well as the 
average number of topplings in a cluster of a specified length I, 

<n> l ~ l Dn (9) 

The scaling exponents are given in Table 1, and as we show in Section 5, they satisfy 
various scaling relations to within numerical error. 

It appears that the exponents are independent of the concentration of defect bonds 
c, suggesting universal criticality. Similar robustness is exhibited by the system for 
variations in A, provided that A is strictly smaller than one: namely, all curves for A<1 
coincide with the open circle curves in Fig. 3. 

4 Nonuniversal criticality: random-slope model (C) 

Model C combines the critical height model with a critical slope instability criterion at 
defect sites. At these sites, which are randomly distributed with relative concentration 
c, large columns of grains will form; these relax according to an alternative set of rules: 
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if at least one of the two local slopes downstream of site exceeds a critical value a c 
then toppling occurs towards these neighbors. 

If a k (i,j) = h(i,j)-h(i + l,j k )>a c , (10) 

then 

h(i,j)^h(i,j)-l] h(i + l,j h )^h(i + l,j h ) + l, (11) 

where the index k = ± stands for right (+) or left (— ) forward neighboring site. These 
rules are applied repeatedly until all slopes become subcritical. In order that topplings 
according to the critical height rule, where we set h c = 2, do not affect those topplings 
that proceed through critical slope dynamics, we set a c = 8 h c . After some transient 
time local slopes cr(i,j) ~ o c — 1 are formed at randomly distributed defect sites, while 
the rest of the system topples according to the critical height rule. Due to the nonlocal 
character of the critical slope rule in Eqs. ( [TU| ) - (|TT|), topplings at defect sites may affect 
stability at their upstream neighbors too, and these sites might then topple in the next 
time step. In this way an internal time scale is introduced, although the macroscopic 
transport direction remains only one way, as in models A and B. 

Results of the numerical simulations for the distributions of duration Pit > t ), size 
D(s > So), and length P(l) are shown in Fig. f| for three concentrations of defect sites 
c = 0.05, 0.2, and 0.6. The exponents for this model, 9, r, and a defined via Pit > 
to) ~ t~ e , D(s > so) ~ s _r , and P(l) ~ respectively, have a (weak) concentration 

dependence. For instance, 9 = 0.516 ± 0.002 for c = 0.05 increases to 0.581 ± 0.001 for 
c = 0.6. Similarly, a = 0.499 ±0.012 at c = 0.05 changes to a = 0.563 ±0.024 at c = 0.6, 
and r = 0.354 ± 0.001 to r = 0.427 ± 0.001 in the same region. For c > 0.7 the curves 
exhibit a finite curvature due to the multiple topplings. (The exponents given in Table 
1 are for c = 0.2.) 

The average duration of avalanches of selected length /, <£>/ in Fig. |]d, exhibits scale 
invariance, in agreement with the scaling properties of the distributions, namely 

<t> t ~l z , (12) 

where z is the dynamic exponent. We find from Fig. ^d that z = 1.013 ± 0.001 for 
c = 0.6. The dynamic exponent differs from 1 as a consequence of the nontrivial time 
scale in this model, although the values of the other exponents are close to (but, as our 
numerical results suggest, distinct from) the ones in the absence of defects. However, 
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nonuniversal properties such as outflow current are different in the two cases: there are 
no apparent finite-size scaling effects in the presence of defects. In the limit c — > 1, all 
sites become subject to the critical slope rule, Eqs. (|T0| ) - (|iT|), due to which a finite 
slope is formed. SOC is lost, since every avalanche is of infinite duration: the dynamics 
is dominated by single grain (one in/one out) events. 



8 



5 Nonequilibrium phase transitions 



In the preceding sections we have shown that sandpile automata are able to self-organize 
into a critical state in the presence of frozen-in random defects, provided that the dy- 
namics conserves the number of grains at each time step. The sets of critical exponents 
for all three models are summarized in Table 1. For model A without defects the expo- 
nents are known exactly @. Model C exponents are for concentration c = 0.2 of defect 
sites. (Note that in models A and B time scale is measured in units of length, therefore 
the corresponding exponents are equivalent, i.e., a = 9.) Also given is the mass-to-scale 
ratio D\\ defined with respect to the length parallel to the transport direction 

<s>i ~ l D " , (13) 

where <s>/ is the time-averaged size of clusters of selected length I. In Fig. ^ <s>i 
is obtained from separate numerical simulations for models A, B, and C. The exponent 
D Tl ( cf. Eq. (H) ) is also given in Table 1 . The numerical values of the exponents are in 
reasonable agreement with the following scaling relation which can be shown to hold in 
all directed models: 

6 z = D\\ r = D n r n = a. (14) 
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Figure 5: Average size of the relaxation clusters <s>z of selected length /, measured 
parallel to the transport direction, for models A, B, and C, as indicated. 
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Table 1: Critical exponents and mass-to-scale ratio in 2D directed models with defects 



Model 


D\\ 


9 


T 


a 


z 









D n 


Remark 


A 


3/2 


1/2 


1/3 


1/2 


1 


1/2 


1 


1/3 


3/2 


universal 


B 


1.62 


1.04 


0.65 


1.04 


1 


none 


1.006 


0.466 


1.997 


universal 


C 


1.49 


0.57 


0.38 


0.54 


1.002 


none 


0.978 


0.366 


1.54 


nonuniversal 



An order parameter which is appropriate for the defect-mediated phase transitions which 
occur in all three models can be defined as q(c) = 1 — < J(c)>/< J(0)>, where <J(c)> the 
outflow current, which is the total number of particles that flow over the lower boundary 
of the system in the presence of defects of concentration c. <> denotes an average over 
the total number of Monte Carlo steps. 

2.0 i 1 1 1 1 1 1 1 1 . 
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Figure 6: Order parameter q(c) vs. concentration of defects c for models A, B, and C, 
as indicated. 



In the models with conservative dynamics, i.e., model B with A = 1 and model C, 
the average flux out of the system is balanced by the incoming current, thus leading 
to the vanishing of the order parameter in the self-organized critical state, as shown 
in Fig. 6. In model A, however, < J(c)>/< J(0)> decreases due to nonconservation of 
number of particles at defect sites, leading to the appearance of finite order parameter q 



10 



with increasing concentration of defect sites. For c > c* = 0.295 = 1 — pd, the directed 
percolation threshold, the states are such that only finite avalanches occur, while for c<c* 
there are relaxation clusters of all sizes. The slope of the q(c) curve at c = c* is close to 
zero. In this way, model A is subcritical for finite concentration c, which appears as a 
control parameter tuning the coherence length of the self-organized state. In the limit 
c — > the coherence length diverges, as discussed in Sec. 2. Our numerical values indicate 
that the order-parameter curve approaches the vertical axis with a large but finite slope. 
We find the same phenomenon in model B with nonconserving dynamics (A<1), where 
transfer is incomplete at each negative bond. (Model B is symmetric with respect to 
transformation c —>■ 1 — c, reflecting the symmetry between positive and negative bonds, 
as can be seen in Fig. 6.) 

There is another type of defect-mediated phase transition in model C in the limit c — > 1. 
At this point SOC is lost in favor of a state with finite net slope. This nonequilibrium 
phase transition requires a more detailed study. In the related case of annealed defects- 
when the probability of toppling p varies at each time step but is the same for all sites 
in the system- a collective phase transition appears at p c = 0.293 [ |TT|| . 



6 Summary 

The presence of frozen-in defects in two-dimensional directed sandpile automata leads to 
various new phenomena. In this paper, we have introduced and studied three variations 
of the directed abelian sandpile automaton on the square lattice in order to explore 
this problem. We find that for site disorder, if the dynamics is nonconservative at 
defects, the system is driven into a subcritical state with a finite correlation length, which 
depends on the concentration c of defects. With any concentration of bond disorder, 
the correlation length remains finite and independent of the degree of nonconservation, 
provided that there is some loss of conservation, A<1. (Other models of nonconservative 
cellular automata have been studied recently |12|], and are seen to have robust SOC 



behavior. Here we have lack of conservation occurring solely due to the presence of 
defects.) Furthermore, if the dynamics is conservative, the automaton self-organizes into 
a critical state with universal scaling properties. For a third type, the case of random 
slope defects, which correspond most closely to the random field version of analogous 
spin problems, we observe some indications of nonuniversal, concentration-dependent 
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scaling exponents. Scaling relations between the exponents are fulfilled exactly for the 
first model (with c = 0), and within numerical error for the latter two models for each 
concentration of disorder. Finally, varying the concentration of defects appears to be a 
mechanism for continuously tuning the local rules of relaxation, which may eventually 
lead to a phase transition between metastable states with different properties. 
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